///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file molpro_interface.cpp
///
/// \brief   contains definitions of Molpro access functions
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#include <src/pcps.h>
#include <src/input.h>
#include <src/hamiltonian.h>
#include <src/blas_lapack.h>
#include <src/molpro_interface.h>

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   reads the number of orbitals from the Molpro FCI dump file
///
/// \param[in]   userinp  user input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////

int pcps::get_molpro_norb(const pcps::Input & userinp) {

  // get molpro dump file name
  const std::string filename = userinp.workdir() + "FCIDUMP";

  // open molpro dump file
  std::ifstream ifs(filename.c_str());

  // check that file is open
  if (!ifs.is_open())
    throw pcps::Exception("failed to open Molpro dump file %s") % filename;

  // skip to the number of orbitals
  std::string skip_str;
  ifs >> skip_str >> skip_str;

  // read the number of orbitals
  int n;
  ifs >> n;

  // close the file
  ifs.close();

  // return the number of orbitals
  return n;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   reads the one and two electron integrals from the Molpro FCI dump file
///
/// \param[in]   userinp  user input options
/// \param[out]  enuc     the nuclear repulsion energy
/// \param[out]  oei      the one electron integrals, size = na*na
/// \param[out]  tei      the two electron integrals, size = na*(na+1) * ( na*(na+1) + 2 ) / 8
///
///////////////////////////////////////////////////////////////////////////////////////////////////

void pcps::get_molpro_ints(const pcps::Input & userinp, double & enuc, double * const oei, double * const tei) {

  // initialize the nuclear repulsion energy
  enuc = 0.0;

  // get numbers of orbitals
  const int no = pcps::get_molpro_norb(userinp); // total
  const int nc = userinp.n_cor_orb();            // core
  const int na = userinp.n_act_orb();            // active

  // check that the numbers of orbitals are sane
  if ( nc + na > no )
    throw pcps::Exception("too many core and/or active orbitals.  nc = %i, na = %i, no = %i") % nc % na % no;

  // check that the number of active orbitals matches the number of lattice sites
  if ( userinp.nsites() != 2*na )
    throw pcps::Exception("number of active orbitals must equal lattice sites:  na = %i, nsites = %i") % na % (userinp.nsites()/2);

  // create temporary vectors into which we will load the integrals
  std::vector<double> temp_int1(no*no, 0.0);
  std::vector<double> temp_int2( no*(no+1)/2 * ( no*(no+1)/2 + 1 ) / 2, 0.0);

  // get molpro dump file name
  const std::string filename = userinp.workdir() + "FCIDUMP";

  // open molpro dump file
  std::ifstream ifs(filename.c_str());

  // check that file is open
  if (!ifs.is_open())
    throw pcps::Exception("failed to open Molpro dump file %s") % filename;

  // skip header
  std::string skip_str;
  while (skip_str != "&END")
    ifs >> skip_str;

  // read integrals
  while (true) {

    // create variables to hold one line's integral information
    double val;
    int i, j, k, l;

    // read a line
    ifs >> val >> i >> j >> k >> l;
    i--; j--; k--; l--;

    // check if we have reached the end of the file
    if (ifs.eof()) break;

    // nuclear energy
    if ( i == -1 && j == -1 && k == -1 && l == -1 )
      enuc += val;

    // one electron integrals
    else if ( k == -1 && l == -1 )
      temp_int1[no*i+j] = temp_int1[no*j+i] = val;

    // two electron integrals
    else
      temp_int2[pcps::Hamiltonian<double>::tei_compound_index(i, j, k, l)] = val;

  }

  // close the file
  ifs.close();

  // extract active space one electron integrals
  for (int p = 0; p < na; p++)
  for (int q = 0; q < na; q++) {
    oei[na*p+q] = temp_int1[no*(p+nc)+(q+nc)];
    for (int i = 0; i < nc; i++)
      oei[na*p+q] +=   2.0 * temp_int2[pcps::Hamiltonian<double>::tei_compound_index(i, i, p+nc, q+nc)]
                     - 1.0 * temp_int2[pcps::Hamiltonian<double>::tei_compound_index(i, p+nc, i, q+nc)];
  }

  // extract active space two electron integrals
  for (int p = 0; p < na; p++)
  for (int q = 0; q <= p; q++)
  for (int r = 0; r < na; r++)
  for (int s = 0; s <= r; s++) {
    const int tmp_cmpd = pcps::Hamiltonian<double>::tei_compound_index(p+nc, q+nc, r+nc, s+nc);
    const int act_cmpd = pcps::Hamiltonian<double>::tei_compound_index(p, q, r, s);
    tei[act_cmpd] = temp_int2[tmp_cmpd];
  }

  // add core orbital contribution to the nuclear energy (nuclear energy is used to also hold core energy)
  for (int i = 0; i < nc; i++)
    enuc += 2.0 * temp_int1[no*i+i] + temp_int2[pcps::Hamiltonian<double>::tei_compound_index(i, i, i, i)];
  for (int i = 0; i < nc; i++)
  for (int j = 0; j <  i; j++)
    enuc +=   4.0 * temp_int2[pcps::Hamiltonian<double>::tei_compound_index(i, i, j, j)]
            - 2.0 * temp_int2[pcps::Hamiltonian<double>::tei_compound_index(i, j, i, j)];

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   throws an error because Molpro doesn't use complex integrals
///
/// \param[in]   userinp  user input options
/// \param[out]  enuc     the nuclear repulsion energy
/// \param[out]  oei      the one electron integrals, size = na*na
/// \param[out]  tei      the two electron integrals, size = na*(na+1) * ( na*(na+1) + 2 ) / 8
///
///////////////////////////////////////////////////////////////////////////////////////////////////


void pcps::get_molpro_ints(const pcps::Input & userinp, double & enuc, std::complex<double> * const oei,
                                                                     std::complex<double> * const tei) {

  throw pcps::Exception("Cannot read complex integrals from Molpro.  Please use scalar_type = real");

}
